/*
 * OptimizationObject.cpp
 *
 *      Author: Jerome RD Soine and Christoph A Brand
 *      Institution: Schwarz goup, Institute for Theoretical Physics, Heidelberg University, Germany
 */

#include "OptimizationObject.h"
#include "include/core/comparator/L2Comparator.h"
#include "include/core/worker/WorkerManager.h"
#include "include/base/cell_models/ACM/ActiveCableModel.h"


/*
 * Constructor for use with parameter file.
 */
OptimizationObject::OptimizationObject(string parameter_file)
{
	cutting_point1 = cutting_point2 = vector<double>(3,0);
	parameter_handler = new ParameterHandler();
	declare_parameter_handler();
	parameter_handler->read_input(parameter_file.c_str());
	define_parameters_from_handler();

	mesh_point1 = new Point<3> (-boundary_extension,-boundary_extension,-substrate_depth);
	mesh_point2 = new Point<3> (image_width + boundary_extension,image_height + boundary_extension,0);

	print_parameters();

}

OptimizationObject::OptimizationObject(string pf, string cmf, string df, double y, double p, Point<3> ms, unsigned int degree_of_refinement)
{
	project_folder = pf;
	cell_mesh = cmf;
	disp_file = df;
	young_mod = y;
	poisson_ratio = p;
	mesh_size = ms;
	deg_of_ref = degree_of_refinement;
}

void OptimizationObject::initializeWorkerManager(string acm_filename, string acm_output_folder, vector<double> parameter){
	cout << "\tGenerating ACM\n";
	ActiveCableModel * acm = new ActiveCableModel(acm_filename,acm_output_folder,50.0,parameter);
	cout << "\tGenerating FEM\n";
	FEMCalc * fem = new FEMCalc(get_young_mod(),get_poisson_ratio(),get_mesh_point1(),
			get_mesh_point2(),get_degree_of_refinement(),acm->getForceTransmissionShapes());
	cout << "\tCalling WorkerManager::initialize\n";
	WorkerManager::initialize(acm, fem, new L2Comparator(new DataContainer(get_disp_file())),8);
}

OptimizationObject::~OptimizationObject()
{
	delete parameter_handler;
	delete mesh_point1;
	delete mesh_point2;
}

string OptimizationObject::get_project_folder()
{
	return project_folder;
}

string OptimizationObject::get_cell_mesh_file()
{
	return cell_mesh;
}

string OptimizationObject::get_disp_file()
{
	return disp_file;
}

double OptimizationObject::get_young_mod()
{
	return young_mod;
}

double OptimizationObject::get_poisson_ratio()
{
	return poisson_ratio;
}

Point<3> OptimizationObject::get_mesh_size()
{
	return mesh_size;
}

unsigned int OptimizationObject::get_n_parameters()
{
	return 6;
}

unsigned int OptimizationObject::get_degree_of_refinement()
{
	return deg_of_ref;
}

int OptimizationObject::get_optimization_method()
{
	return optimization_method;
}

double OptimizationObject::get_initial_temperature()
{
	return initial_annnealing_temp;
}

double OptimizationObject::get_simplex_step_width()
{
	return simplex_step_width;
}

double OptimizationObject::get_annealing_parameter()
{
	return annealing_parameter;
}

void OptimizationObject::define_parameters_from_handler()
{
	project_folder = parameter_handler->get("project folder");
	cell_mesh = parameter_handler->get("cell mesh file");
	disp_file = parameter_handler->get("displacement file");
	shape_file = parameter_handler->get("shape file");
	young_mod = parameter_handler->get_double("youngs modulus");
	poisson_ratio = parameter_handler->get_double("poisson ratio");
	reg_param = parameter_handler->get_double("regularization parameter");
	image_width = parameter_handler->get_integer("image width");
	image_height = parameter_handler->get_integer("image height");
	substrate_depth = parameter_handler->get_integer("substrate depth");
	boundary_extension = parameter_handler->get_integer("boundary extension");
	preshifting_enabled = parameter_handler->get_bool("pre-shifting enabled");
	preconditioning_enabled = parameter_handler->get_bool("preconditioning");
	simulated_data = parameter_handler->get_bool("simulated data");
	n_optimization_steps = parameter_handler->get_integer("n optimization steps");
	string optimization_method_string = parameter_handler->get("optimization method");
	if(optimization_method_string == "conjugated gradients") optimization_method = METHOD_CONJUGATED_GRADIENT;
	else if(optimization_method_string == "conjugated gradients nt first") optimization_method = METHOD_CONJUGATED_GRADIENT_NT_FIRST;
	else if(optimization_method_string == "conjugated gradients nt last") optimization_method = METHOD_CONJUGATED_GRADIENT_NT_LAST;
	else if(optimization_method_string == "simplex") optimization_method = METHOD_SIMPLEX;
	else if(optimization_method_string == "simplex annealing") optimization_method = METHOD_SIMPLEX_ANNEALING;
	else if(optimization_method_string == "fixed network tension cg") optimization_method = METHOD_FIXED_NT_CG;
	else if(optimization_method_string == "sample around network tension") optimization_method = METHOD_SAMPLE_AROUND_NT_CG;
	else optimization_method = METHOD_STEEPEST_DESCENT;
	deg_of_ref = parameter_handler->get_integer("n of mesh refinements");
	n_threads = (unsigned int) parameter_handler->get_integer("n threads");
	tension_input = parameter_handler->get_bool("tension input");
	tension_file = parameter_handler->get("tension file");
	n_sf_steps_per_network_tension_adjustment = parameter_handler->get_integer("n sf steps per network tension adjustment");
	generate_results = parameter_handler->get_bool("generate results");
	initial_annnealing_temp = parameter_handler->get_double("initial temperature");
	annealing_parameter = parameter_handler->get_double("annealing parameter");
    simplex_step_width = parameter_handler->get_double("simplex step width");
	known_network_tension = parameter_handler->get_double("known network tension");
	decrease_network_tension = parameter_handler->get_bool("decrease network tension");
	cutting_point1[0] = parameter_handler->get_double("cutting p1x");
	cutting_point1[1] = parameter_handler->get_double("cutting p1y");
	cutting_point1[2] = parameter_handler->get_double("cutting p1z");
	cutting_point2[0] = parameter_handler->get_double("cutting p2x");
	cutting_point2[1] = parameter_handler->get_double("cutting p2y");
	cutting_point2[2] = parameter_handler->get_double("cutting p2z");
}

void OptimizationObject::declare_parameter_handler()
{
    parameter_handler->declare_entry ("cell mesh file",
        "cell.vtk", Patterns::Anything(), "Name of cell mesh file generated with IJ plugin.");

    parameter_handler->declare_entry ("project folder",
        "project", Patterns::Anything(), "Defines folder for saving results.");

    parameter_handler->declare_entry ("displacement file",
        "displacements.txt", Patterns::Anything(), "Defines the file that contains the measured displacements.");

    parameter_handler->declare_entry ("shape file",
        "shapes.txt", Patterns::Anything(), "Defines the file that contains shapes of focal adhesions.");

    parameter_handler->declare_entry ("youngs modulus",
        "1030.4", Patterns::Double(), "Youngs modulus of the substrate nN/px^2");

    parameter_handler->declare_entry ("poisson ratio",
        "0.49", Patterns::Double(), "poisson ratio of the material.");

    parameter_handler->declare_entry ("regularization parameter",
            "0.0", Patterns::Double(), "Defines regularization parameter ");

    parameter_handler->declare_entry ("image width",
         "1200", Patterns::Integer(), "Width of the image in px.");

    parameter_handler->declare_entry ("image height",
         "1200", Patterns::Integer(), "height of the image in px.");

    parameter_handler->declare_entry ("substrate depth",
         "300", Patterns::Integer(), "substrate depth in px.");

    parameter_handler->declare_entry ("boundary extension",
         "300", Patterns::Integer(), "Extension of the boundary in order to avoid effects from zero stress BCs.");

    parameter_handler->declare_entry ("n of mesh refinements",
         "6", Patterns::Integer(), "Number of mesh refinements outgoing from a single cube.");

    parameter_handler->declare_entry ("pre-shifting enabled",
         "true", Patterns::Bool(), "Previous shifting of the mesh based on the displacements information.");

    parameter_handler->declare_entry ("preconditioning",
         "true", Patterns::Bool(), "Preconditioning of the starting vector based line minimization along type directions.");

    parameter_handler->declare_entry ("simulated data",
         "false", Patterns::Bool(), "If true, data is not taken from displacement file but is simulated.");

    parameter_handler->declare_entry ("tension input",
         "false", Patterns::Bool(), "If true, tension file is evaluated as starting point for optimization.");

    parameter_handler->declare_entry ("tension file",
        "", Patterns::Anything(), "Name of tension file.");

    parameter_handler->declare_entry ("n optimization steps",
         "10", Patterns::Integer(), "Number of optimization steps.");

    parameter_handler->declare_entry ("optimization method",
         "steepest descent", Patterns::Anything(), "Name of the optimization method.");

    parameter_handler->declare_entry ("n threads",
         "8", Patterns::Integer(), "Number of threads.");

    parameter_handler->declare_entry ("n sf steps per network tension adjustment",
         "10", Patterns::Integer(), "Number of stress fiber optimizations per network tension adjustment step.");

    parameter_handler->declare_entry ("generate results",
         "false", Patterns::Bool(), "Only generates results. Needs tension input file.");

    parameter_handler->declare_entry ("initial temperature",
            "100.0", Patterns::Double(), "Defines intial temperature for simplex annealing method");

    parameter_handler->declare_entry ("annealing parameter",
                "1.0", Patterns::Double(), "Defines the annealing scheme");

    parameter_handler->declare_entry ("simplex step width",
                    "0.1", Patterns::Double(), "Defines the size of the simplex moves");

    parameter_handler->declare_entry ("known network tension",
        "0.0", Patterns::Double(), "pregiven network tension in combination with mode *fixed network tension cg*");

    parameter_handler->declare_entry ("decrease network tension",
         "true", Patterns::Bool(), "for network tension sampling: true = decrease nt, false = increase nt");

    parameter_handler->declare_entry ("cutting p1x",
         "0.0", Patterns::Double(), "x coordinate of first point that defines a cutting line");

    parameter_handler->declare_entry ("cutting p1y",
         "0.0", Patterns::Double(), "y coordinate of first point that defines a cutting line");

    parameter_handler->declare_entry ("cutting p1z",
         "0.0", Patterns::Double(), "z coordinate of first point that defines a cutting line");

    parameter_handler->declare_entry ("cutting p2x",
         "0.0", Patterns::Double(), "x coordinate of first point that defines a cutting line");

    parameter_handler->declare_entry ("cutting p2y",
         "0.0", Patterns::Double(), "y coordinate of first point that defines a cutting line");

    parameter_handler->declare_entry ("cutting p2z",
         "0.0", Patterns::Double(), "z coordinate of first point that defines a cutting line");
}

void OptimizationObject::print_parameters()
{
	cout << "Parameter Settings:" << "\n";
	cout << "------------------------------------------- \n";
	cout << "Project Folder     = " << get_project_folder() << "\n";
	cout << "Cell mesh file     = " << get_cell_mesh_file() << "\n";
	cout << "Displacement file  = " << get_disp_file() << "\n";
	cout << "Shape file 		= " << get_shape_file() << "\n";
	cout << "Youngs modulus     = " << get_young_mod() << "\n";
	cout << "Poisson ratio      = " << get_poisson_ratio() << "\n";
	Point<3> p = get_mesh_point1();
	cout << "Mesh Point 1       = (" << p(0) << "," << p(1) << "," << p(2) << ")\n";
	p = get_mesh_point2();
	cout << "Mesh Point 2       = (" << p(0) << "," << p(1) << "," << p(2) << ")\n";
	cout << "N of Mesh Refs.    = " << get_degree_of_refinement() << "\n";
	cout << "N of threads       = " << get_n_threads() << "\n";
	cout << "Pre-shifting       = " << is_preshifting_enabled() << "\n";
	cout << "Preconditioning    = " << is_preconditioning_enabled() << "\n";
	cout << "simulated data     = " << is_simulated_data_enabled() << "\n";
	cout << "N of Opti steps    = " << get_n_optimization_steps() << "\n";
	cout << "Method	            = " << get_optimization_method() << "\n";
	cout << "Use tension input  = " << is_tension_input_enabled() << "\n";
	cout << "Tension file       = " << get_tension_file() << "\n";
	cout << "n sf steps per nto = " << get_n_sf_steps_per_network_tension_adjustment() << "\n";
	cout << "Reg. parameter     = " << get_reg_param() << "\n";
	cout << "Initial temperature	= " << get_initial_temperature() << "\n";
	cout << "Annealing parameter 	= " << get_annealing_parameter() << "\n";
	cout << "simplex step width	= " << get_simplex_step_width() << "\n";
	cout << "known network tension	= " << get_known_network_tension() << "\n";
	cout << "decrease network tension	= " << is_decrease_network_tension() << "\n";
	cout << "cutting point 1	= (" << cutting_point1[0] << "," << cutting_point1[1] << "," << cutting_point1[2] << ")\n";
	cout << "cutting point 2	= (" << cutting_point2[0] << "," << cutting_point2[1] << "," << cutting_point2[2] << ")\n";
	cout << "------------------------------------------- \n \n \n";

}

void OptimizationObject::save_parameters(string filename)
{
	ofstream data (filename.c_str());

	data << "Parameter Settings:" << "\n";
	data << "------------------------------------------- \n";
	data << "Project Folder     = " << get_project_folder() << "\n";
	data << "Cell mesh file     = " << get_cell_mesh_file() << "\n";
	data << "Displacement file  = " << get_disp_file() << "\n";
	data << "Shape file 		= " << get_shape_file() << "\n";
	data << "Youngs modulus     = " << get_young_mod() << "\n";
	data << "Poisson ratio      = " << get_poisson_ratio() << "\n";
	Point<3> p = get_mesh_point1();
	data << "Mesh Point 1       = (" << p(0) << "," << p(1) << "," << p(2) << ")\n";
	p = get_mesh_point2();
	data << "Mesh Point 2       = (" << p(0) << "," << p(1) << "," << p(2) << ")\n";
	data << "N of Mesh Refs.    = " << get_degree_of_refinement() << "\n";
	data << "N of threads       = " << get_n_threads() << "\n";
	data << "Pre-shifting       = " << is_preshifting_enabled() << "\n";
	data << "Preconditioning    = " << is_preconditioning_enabled() << "\n";
	data << "simulated data     = " << is_preconditioning_enabled() << "\n";
	data << "N of Opti steps    = " << get_n_optimization_steps() << "\n";
	data << "Method	            = " << get_optimization_method() << "\n";
	data << "Use tension input  = " << is_tension_input_enabled() << "\n";
	data << "Tension file       = " << get_tension_file() << "\n";
	data << "n sf steps per nto = " << get_n_sf_steps_per_network_tension_adjustment() << "\n";
	data << "known network tension	= " << get_known_network_tension() << "\n";
	data << "decrease network tension	= " << is_decrease_network_tension() << "\n";
	data << "cutting point 1	= (" << cutting_point1[0] << "," << cutting_point1[1] << "," << cutting_point1[2] << ")\n";
	data << "cutting point 2	= (" << cutting_point2[0] << "," << cutting_point2[1] << "," << cutting_point2[2] << ")\n";
	data << "------------------------------------------- \n \n \n";

	data.close();
}
